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Abstract We study collective dynamics of complex networks of stochastic excitable 
elements, active rotators. In the thermodynamic limit of infinite number of elements, 
we develop a mean-field theory for the network and then use a Gaussian approximation 
to obtain a closed sequence of deterministic integro-differential equations. These equa- 
tions govern the order parameters of the network. We find that a homogeneous decrease 
in the number of elements connections merely shifts the bifurcation thresholds without 
producing qualitative changes in the network dynamics. In contrast, heterogeneity in 
the number of connections leads to bifurcations in the excitable regime. In particular 
we show that a critical value of noise intensity for the saddle-node bifurcation decreases 
with growing connectivity variance. The threshold values of noise intensity for the on- 
set of global oscillations (Hopf bifurcation) show a non-monotone dependency on the 
structural heterogeneity, displaying a minimum at moderate connectivity variances. 

1 Introduction 

Excitable elements are ubiquitous in physical, geophysical, chemical and biological systems. Exam- 
ples range from chemical reactions, neuronal systems, cardiovascular tissues and climate dynamics 
to laser devices, see Ref. |1I2) and references therein. An excitable system possesses a stable equilib- 
rium, but, once excited by sufficiently strong perturbations, displays large non-monotonic excursions 
before returning to the state of rest. 

Networks of excitable elements are omnipresent in nature, and the human brain, as a gigantic 
network of neurons, appears to be paradigmatic in this respect (cf. ^). In networks of nonlinear os- 
cillators, one typically encounters the collective phenomenon of synchronization: different oscillators 
adjust their rhythms due to coupling 4 5 6J. Recently much effort has been invested into under- 
standing the mechanisms which stand behind synchronization in complex networks of nonlinear 
oscillators and excitable systems, see Ref. |7I8| and references therein. 

A crucial step in characterizing limit cycle oscillators was done by Kuramoto 'W : for a broad class 
of systems, he reduced description to dynamics of a single scalar variable, the uniformly rotating 
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phase. Since the origins of this approach can be traced back to the biologically motivated Winfree 
model |10j , it appears to be a good candidate for providing fundamental insights into the function- 
ing of coupled biological oscillators. In general, coupled limit-cycle oscillators can be approximated 
by the Kuramoto model as long as the interaction between them stays weak, so that each oscil- 
lator experiences only a small phase perturbation. This circumstance explains significance of the 
Kuramoto approach and its modifications in the context of computational neuroscience |11I12| . Rel- 
ative simplicity of the Kuramoto model turns it into a convenient tool for analytical explorations. 
In the frame of this model, characteristics of different transitions between desynchronized and syn- 
chronous states in large ensembles of oscillators can be obtained in closed form. This explains 
continuing interest in the Kuramoto model among physicists, neuroscientists and mathematicians 

mm- 

Phase oscillators in their canonical form are not directly appHcable to studies of coupled excitable 
elements. In particular a Kuramoto oscillator possesses translational invariance with respect to 
phase shifts: it rotates uniformly, and all phase values are dynamically equivalent, whereas for an 
excitable element a certain event (e.g. spike) is singled out. An "active rotator" ~ a modification 
of the Kuramoto model which takes into account these properties of excitable elements - was 
introduced by Shinomoto and Kuramoto in [15j . Since then, the analysis was refined and certain 
additional aspects were investigated in [16117118119120] . Most of the research considered the case 
of uniform global coupling in the ensemble; only in a few works, active rotators were studied on 
a network with complex structure, see e.g. |21] . Here we study dynamics of networks of randomly 
connected stochastic active rotators and put emphasis on how temporal fluctuations and the network 
structure influence the dynamics of the network as a whole. 

The paper is organized as follows. In section [2 we introduce the basic model and briefly discuss 
its properties. In section [3] the mean- field approximation is developed: the Fokker-Planck equation 
for the time-dependent probability distribution of the phases is derived and transformed by means 
of the Fourier expansion into the infinite hierarchy of ordinary differential equations. In section |3] 
this hierarchy is brought to a closed form with the aid of a Gaussian approximation. Our approach 
is exemplified for regular and binary random networks in section [S] Numerical simulation of the 
binary network corroborates the conclusions drawn from the analysis of the mean-field model. 



2 Model 

In the population of noise-driven active rotators, introduced by Shinomoto and Kuramoto [15j . 
dynamics of individual phases (/'i(t) is governed by equations 

N 

Mt) = 1 - a Sin (Mt)) + ^Y. - '^»(*)) + ^»(*)' 

i=i 

where the units are indexed by i = 1, . . . , iV. Parameter a determines the excitation threshold of 
each separate element, whereas k characterizes the coupling strength. Here, we restrict ourselves 
to the case where these parameters are identical for all elements of the ensemble. The network is 
assumed to be undirected and unweighted, hence the elements of the symmetric adjacency matrix 
A adopt only two values: Aij = 1, if the units i and j are coupled, and Aij = otherwise. The 
number of connections of a node, given by 



N 

h = ^ Aij , 



(2) 
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Figure 1. (Color online) Effect of the mean- field approximation on an undirected network of four nodes. 
Left: original complex network. Right: the approximate network. Thickness of the edges is approximately 
proportional to their weight. 



is called the degree of this node. In absence of isolated nodes, ki > 0. Complex topologies of real- 
world networks can be encoded into the adjacency matrix, and decoded by counting all the degrees, 
which gives rise to a degree distribution P{k). Throughout this paper we discuss random networks. 

The noise terms (t) in ([T]) are statistically independent across the network nodes and are zero 
mean Gaussian white sources, {£,i{t)) = 0, = '2D5ijS{t — t'), where the angular 

brackets denote averages over different realizations of the noise and D is the noise intensity. 

For an isolated element in the absence of noise, (pi = 1 — asin excitable behavior is apparent. 
For \a\ > 1 the stable equilibrium is at = arcsin(a~^), and the system needs a sufficiently strong 
perturbation in order to make an excursion around the circle. Noise plays this role driving the system 
to escape from the equilibrium state <j>°° . An escape event corresponds to the release of a single 
spike [1512211] . For |a| < 1 the element shows oscillatory behavior with frequency Vl — a^- Notably, 
evolution of the phase (fi is not uniform in this case: it is at the slowest near (jj — tt/2 and at the 
fastest near </> = 37r/2. 



3 Mean-field theory 



A mean-field description can be constructed by assuming that all nodes in the network are statisti- 
cally identical, so that they are fully characterized by their individual degrees [23] . The derivation 
of the mean-field description can be done by coarse-graining the network structure, namely by re- 
placing the original network by a fully connected one with random coupling weights that mimic 
the actual connectivity (see Fig. [T] for illustration) . Considering an uncorrelated random network 
[24|23j25] and requiring the conservation of the individual degrees, ki = X^jLi ^iji i — 1, . . . , N , 
the elements of the approximate adjacency matrix read 

^ij = h jv^ ■ (3) 
Ei=i ki 

Inserting this into Eq. ([T|) yields a weighted mean-field approximation and effectively provides a 
one-oscillator description. 

The system of N coupled Langevin equations ([T|) can be characterized in terms of the transition 
probability density Vn {(p, t\<t>^, fe) . Here, the vector <f) = {cpi is built from the phases 
of N rotators at time t whereas (jP = {4>\t . . consists of their values at initial time 

Components of the time-independent vector k = (fci, . . . , kjq) are the degrees of all nodes in the 
network. We assume that the initial phases 0^ are identically and independently distributed. The 
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transition probability density obeys a Fokker-Planck equation (FPE) which describes the evolution 
of the ensemble from time t° to time t > t°. 

Integration of the FPE over the iV — 1 phases 02, • ■ ■ , 07v, degrees k2, ■ ■ ■ and corresponding 
initial phases [53] , yields an evolution equation for the one-oscillator probability density Pi (0i , i, | fci ) , 
in which, however, the two-oscillator density P2(0i, 4'2,t, 1^1,^2) is present as well (see Ref. [33] for 
the case of all-to-all coupling with more general coupling functions): 



(4) 



d_ 



J^~^~^T~^ /"d02^sin(02 - 0l) ^2^2 (0i,02,t; 



In the thermodynamic limit N 00 the ratio J^i ki/{N — 1) approaches the average degree (fc'). In 
the framework of the mean-field approximation, correlation of phases between any two oscillators 
can be discarded, 7^2 (0i, 02, fci, fe) ^ "Pi (0i, t; fci) T^i (02, t; ^2) . In this way Eq. ^ becomes 
closed in Pi(0i, t, \ ki), yielding an effective one-oscillator description; the price is, that it becomes 
a nonlinear FPE [27] : 

dt = --Q^[v{(l},t\k)Vi{(l),t\k)] + D — . (5) 

Since all oscillators are statistically identical in this setting, from now on we neglect the indices 
for the individual phases and degrees. The mean increment of the phase per unit time, v{4>,t\k) = 
1 ~ a sin -|- rnkN^^ sin(6' — 0) depends on the density Vi (0, t\k) via the mean-field amplitude r{t) 
and phase 0{t), 

r2iT 



rity""^'^ = / " d0' X] e'*' Vi (0', t\k') k' P [k') , (6) 



where P{k') denotes the degree distribution. We note that Eq. © suggests to introduce order 
parameters for each degree class, namely 

r(Oe^^(*) ^ J_ ^ r,, ityo,, it) ^, p 

^ ' k' 

Let us proceed with the Fourier series expansion of the one-oscillator probability density: 

+00 

Pi(0,t|fc)-— J2 Pn(i|fc)e™* , (8) 

71 — — 00 

with po — 1 p-n — pn- Through the inverse transform of ()15p one can write 

p„(i|fc) = ^ E e"^'' = + ^Snim, (9) 

je{fe} 

where n = 1,2, .. . ,00; Nk stands for the number and {k} for the set of nodes with degree k. For 
every degree k that appears in the network, the Fourier coefficients p„(i|fc) are governed by an 
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infinite chain of coupled complex- valued equations 
Pn{t\k) Rk 



ypn-Mk) ^ p,{t\k')k'p{k') p„+i(i|fc) p^,{t\k')k'p{k') 

+ I (p«-i(i|fc) - Pn+i{t\k)) - (Dn - i)pn{t\k), 



(10) 



with the abbreviation k :— k/{2N). The mean degree (fc') of the network is the first moment of 
the degree distribution. As expected, in the case of all-to-all connectivity the existing results of 
[28j are recovered. Equation pop provides a complete description of our system in the mean-field 
approximation, which is valid as long as A: ^ 1. Thus, for a sparsely connected network the mean- 
field description is likely to fail 29 23 30] . We also note that the above mean field approach is 
restricted to networks with finite second moment of the degree distribution |29|31|23j . 



4 Gaussian approximation 



There seems to be no gain in reduction of complexity from recasting the Langevin dynamics (U)) 
into the mean field description (jlOp . One approach to turn (|10l) into a dynamical system of finite 
order is to truncate the Fourier series. In the case of all-to-all coupling the values of \pn \ decay with 
increase of n sufficiently fast, and the quantitative dynamics is correctly captured by the first two 
dozens of coefiicients so that the rest of the series can be neglected [28 32 . An alternative approach 
lies in finding a framework in which a closure of the infinite set of equations (jlOp can be achieved. 
The so-called Ott-Antonsen ansatz [33; yields impressive results for deterministic ensembles of 
coupled Kuramoto oscillators, reducing exactly the infinite-dimensional set of equations to a low- 
dimensional system. Unfortunately, for the case of noisy systems the direct application of the Ott- 
Antonsen ansatz is impossible, and we are unaware of its appropriate modifications. Instead, we 
use a Gaussian approximation (GA): we assume that in every subset of oscillators with the same 
network degree, the distribution of the phases at every moment of time is Gaussian with the mean 
mkit) and the variance cr^(i). The smaller the noise intensity D, the more accurate is the GA [22) . 
The noise should not dominate the collective dynamics, i.e. D < {k)N~^ . 

In the thermodynamic limit, N ^ oo, the averaging in Eq. ^ in the GA yields 

Cnit\k) = e-"''^'(*)/2 (nrukit)) , s„{t\k) = c-"''^^^/^ sin (nmfe(t)) , (11) 

therefore all c„ and s„ are predetermined by ci and si: C2 — c\ — sf, S2 — 2siCi (sf -I- c^), etc. 
[32] • Thus, the dimension of the reduced system of equations for the Fourier coefficients is twice the 
number of different degrees present in the network. Denoting ci{t\k) by ci and si{t\k) by si, Eqs. 
(fTU]) turn into 

= ^ [{ci{t\k')k'){l -4 + s\) - 2{s^{t\k')k')s^c^{sl + c?)] +^{l-c\ + s\) - Dc, - si, 
^1 = ^ [(5i(i|fc')fc'>(l + Cl - 4) - 2{ci{t\k')k')sici{sl + cD] - asiciisl + c?) - Dsi + ci, 
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Figure 2. (Color online) Bifurcation diagram of a regular network of stochastic rotators in the GA. (a): 
The coupling strength equals k = 1. sni^2'- saddle-node bifurcation lines; H: Hopf bifurcation line; h: 
homoclinic bifurcation, S: oscillation type changes; BT: Bogdanov-Takens bifurcation point, L: homoclinic 
loop attached to the saddle-node, (b): Fragments of the bifurcation diagram for the indicated values of the 
connectivity fraction a. 



for every degree k in the network. Averaging (...) is made over the degree distribution P{k'). 
Equivalently, equations for the mean mk{t) and the variance cr^(i) are 



TOfc =1 — e '^'"^^ cosher^ 



a sin TOfc — — ( (fc'e^'^fc' " sin mk' ) cos mk — {k'e^"h' " cos mk' ) sin rrife ] 

N{k') \ J 

-""'/^sinhCT^ 

nik + , , ( (fc'e"'^fc'/^ sinmfc/) sinmfc + {k'e~'^i" /'^ cos ruk') cos mk) 



(13) 



5 Results: regular and binary random networks 

The integro-difFerential character of Eqs. (|13p makes further general analysis cumbersome. Below 
we discuss explicitly two relatively simple examples of regular and binary random networks. 



5.1 Regular networks 

For regular networks where all nodes have the same degree, Eqs. ([T^ simplify to 



ruk 



= 1- ae "^^^ cosher? sin 



&l = 2D-2e-'''^/^smhal 



a cos TOfc 



+ nkN-^c-'^^^/^ 



(14) 



Here the only difference to the known cumulant equations for a fully connected network [32j is a 
rescaling of the coupling strength: k — > nk/N. Fig. [5] shows a complete bifurcation diagram for 
all-to-all coupling on the parameter plane excitability a and noise intensity D. Two saddle-node 
bifurcations, denoted by srii and sn2, separate parameter regions in which only a single steady 
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Figure 3. (Color online) Critical noise intensity Dc as function of the connectivity fraction a. (a): Small and 
moderate values of a, only the Hopf bifurcations are present, (b): Large values of a with Hopf and saddle- 
node (SN) bifurcations present. Parameter values are indicated in the legends. The coupling strength is 
K. = 1 for both panels. 



state exists from the central wedge-shaped region in which three steady states are present. While 
at sn2 a saddle point and a node are created, sni destroys the saddle and the original steady state. 
On both bifurcation curves one of the Jacobian eigenvalues vanishes. There is a point where the 
second eigenvalue vanishes as well. This codimension-2 Bogdanov-Takens (BT) bifurcation is an 
origin of two further bifurcation lines: the Hopf bifurcation H and the homoclinic bifurcation h. 
While above H, the steady state is stable, below H the stable limit cycle comes into existence. The 
curve h marks the existence of an orbit homoclinic to the saddle point. In the parameter region 
between the curves TJ, h and sn2 the stable steady state coexists with the attracting limit cycle. 
In the region between H and sni_2 two stable steady states exist. At the point L, saddle-node and 
homoclinicity merge together [32] ■ Below the Hopf bifurcation H , small-amplitude oscillations are 
born. With decreasing noise intensity D, their amplitude grows. Below the line S, oscillations turn 
into full rotations: instead of oscillating back and forth between two extremal values, the distribution 
rotates around the whole circle [0, 2tt). 

In the following, we do not discriminate between localized oscillations and full-circle rotations. 
We rather focus on how oscillating patterns are affected by the network structure. Therefore, we 
concentrate on locating the Hopf and the saddle-node bifurcations. When distinguishing between 
oscillating and non-oscillating behavior one has to bear in mind that the two different stationary 
states have different physical nature. Below the saddle-node bifurcation line, the collective dynamics 
is not time-periodic, because the oscillators are mostly at rest without producing any current rrifc ~ 
|17j . On the contrary, above the Hopf bifurcation line, the system is also in a stationary state, but 
with a non-zero current, which grows with further increase of D until saturation. In this state 
individual units are firing, but do not create a globally synchronized oscillation. 

Since a reduction in the number of connections effectively decreases the coupling strength, no 
new qualitative phenomena are expected if all nodes share the same degree. In the following, we 
use the abbreviation a := k/N for the connectivity fraction. 

We detect numerically the dependence on a of Hopf and saddle-node bifurcations in Eqs. 
for different excitation thresholds a with the help of the software package MATCONT 34J . Results 
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are summarized in Fig. [3] showing the bifurcation value of noise intensity (called the critical noise 
intensity) Dc versus a for fixed values of the excitability parameter a. We observe that both the 
Hopf and the saddle-node bifurcations are shifted to smaller noise intensities when the connectivity 
is decreased. Dependence of the Hopf threshold on a is practically linear, in accordance with the 
expansion 

D, = na(--—a^- —a^ + . . .V (15) 
V 2 32 256 y ^ ' 

As one would expect, the oscillatory regions shrink for higher values of threshold a. 



5.2 Binary random networks 



A simple example of a complex heterogeneous network is a binary random network which possesses 
two distinct connectivity degrees. Denoting these as fci and fc2 (ai — ki/N, a2 ~ k2/N), the degree 
distribution is given by 

P{k)=pSkM + i'^-p)3k,k,. (16) 

These networks tend to be disassortative, i.e. nodes with different degrees are preferentially con- 
nected [35]. Such degree-degree correlations could probably be factored in by developing a theory 
along the lines of [35]. It turns out that the disassortativity attains large values, if the variance 
among the connections is large. The maximum variance here is VarMax — 0.25 and we consider 
only values up to Var(Q:) /VarMax = 0.5, where the disassortativity is negligible. Inserting the degree 
distribution (116p into Eqs. (fT3|) gives a system of four coupled first-order differential equations: 

Kaia2(l — p) -a^ /? 

asin(mi)-| — sin(mi — m2)e ^' 

pai + (l-p)a2 



mi =1 - e "'i/^ cosh(cr^ 
&l =2D - 2c-''^^/'^ smh{af) 



1712=1- e "'^^^ cosh((T 



dl =2D - 2c"'"2'/2 sinh(cr^) 



(mi)H — ^ — (aipe '^i/^ + q;2(1 - p)e '^^^^ cos{mi - m2)) 

pai + (1 — p)ct2 ^ ' 

asin(m2)-l -. r — sin(m2 — mi)e ^' 



a cos(m2) + 



Ka2 



pai -I- (1 — p)a2 



pai + (1 - p)a2 
(a2{l -p) 



(17) 



e "'2/2 ^ Q,^pg '^i/'^ cos(m2 — TOi) 



Bifurcation analysis of this system was performed using MATCONT. 

Fig. Ufa) illustrates a binary random network. In (b) we show the order parameters 7'i,2(i) Q 
obtained by simulation of the full dynamics ([1]). More connections appear to increase the similarity 
between the phases in the corresponding population. Fig. \Ms) shows phase portraits of the GA 
Eqs. pT|) near the saddle-node bifurcation. The left panels depict the oscillatory regime. Indeed, 
the set of nodes with more connections is more synchronized both in the case of localized oscillations 
[upper left] and for full rotations [lower left]. This is signalled by smaller variances cr|(t). In the 
right panels of (c) a tiny decrease in the noise intensity brings the whole system into a stationary 
state, where the variables toi,2, cr^ 2 approach a fixed point. 

In this context we never observed that one population of a certain degree is oscillating, while the 
other is not. This might be the result of entanglement of two populations with different connectivity 
degrees. Hence, a visualization such as in Fig. HJa) has to be taken with care. 



Will be inserted by the editor 



9 




5 10 15 20 



Figure 4. (Color online) Dynamics of a binary random network of stochastic rotators, (a): A snapshot 
at i = 100 of a network composed of 50 elements. For clarity, the nodes are distributed on two separated 
rings according to their degrees and phases. Blue nodes have degree fci = 4, red ones have k2 — 20. Green 
links are chosen between different degrees, while the orange lines connect nodes with the same degree. 
Remaining parameters are a — 0.2, D — 0.1, k — 2, p = 0.7. (b): Order parameters ri{t) = |/9i(t|fci) 
and r2{t) — |pi(t|fe2)| obtained from simulation of the full dynamics, Eq. ((l|. The parameters are: N — 
W^, ki = 1000, k2 = 4000, p = 0.7, a = 0.6, D = 0.1, k = 1.5. (c): Phase portraits of the GA Eqs. ^ 
for ai = 0.37, 02 = 0.97, p — 0.95, k = 1 and initial conditions mi_2(t = 0) = o"i,2(^ = 0) = 0.5. Left 
panels depict oscillatory regime, whereas in the right panels a small decrease of noise has the result that 
the system approaches fixed points (right endpoints of the trajectories). 



In order to isolate the effects of network heterogeneity, we fix the average degree, (a) = 0.4. 
The variance of a is tuned by varying three parameters: ai, 02 and p. Since Var(a) is not uniquely 
defined and different parameter sets yield the same variance, we take the path where the a-values 
are the largest, because then the mean-field assumption is best fulfilled. Fig. shows critical noise 
intensity Dc at which a Hopf or saddle-node bifurcations occur in the GA system (fT7| versus the 
connectivity variance Var(a)/VarMax. In the oscillatory regime, a < 1, shown in Fig. ^a), Dc vs 
Var(a)/VarMax shows trends similar to the case of a = 23J, 

^ , „N K (fc'^) n / , , Var(a)\ , „, 

Nonlinear dependence of Dc vs Var(Q!) starts to appear only near the threshold a = 1. In the 
excitable regime a > 1 [Fig. [Sfb)] the oscillatory region grows with the increase of Var(a), if the 
variance is already large. In particular, the Hopf bifurcations shift to larger noise values, while 
the saddle-node bifurcations shift to smaller ones. Furthermore, for the Hopf bifurcation the de- 
pendence Dc vs Var(a) becomes non- monotone with minima at moderate values of the structural 
heterogeneity. These minima become more pronounced for larger values of a. This effect constitutes 
a qualitative difference to the case of oscillatory regime, a < 1, and is induced by the network 
structure. 

In particular, the minima observed for the Hopf bifurcation lines result from adding a small 
fraction of highly connected nodes, so-called hubs. If the nodes with larger degree are in the majority, 
then the critical noise is larger and the minima disappear (not shown here). 
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Figure 5. (Color online) Critical noise intensity versus Var(ci)/VarMax obtained in the GA for the 
binary random network of stochastic rotators. As in Fig. [3] panels (a) and (b) show the oscillatory and the 
excitatory regimes, respectively. The values of excitability parameter a are indicated in the legends along 
with the type of bifurcation: Hopf or saddle- node (SN). Other parameters are VarMax = 0.25, (a) = 0.4, 
K = 1. 

5.3 Numerical simulations of binary random networks 

To test theoretical predictions made with the mean-field theory in Gaussian approximation, we 
performed numerical simulations of random binary networks of stochastic rotators. The goal was 
to check whether the global oscillations can indeed be suppressed by increasing the network het- 
erogeneity from a small to a moderate value. We took networks of = 10** rotators and integrated 
Eq. ([T]) using the Heun method [37] with the time step of 0.5. Smaller time steps down to 0.05 had 
no significant effect on the simulation results. 

We used the Viger-Latapy algorithm to generate the binary random networks |38| . In partic- 
ular, we used the C library igraph [i39j (version 0.6.5), where the aforementioned algorithm is 
implemented. 

Fig. [S] shows simulation results for two different networks. Both networks had the same average 
degree, but distinct variances. Compared with Fig. O the only difference in parameters of Fig. El 
is a larger value of the coupling strength which has no qualitative effect on global dynamics but 
simplifies the detection of the global oscillations. As can be seen from Fig. [51 an increase of the 
network heterogeneity leads to suppression of network's oscillations, in agreement with the theory. 

6 Conclusion 

We have studied the well-known model for excitable dynamics, the so-called active rotators, in the 
context of heterogeneous networks. Besides investigating the effects of noise, we put emphasis on 
the effects that emerge due to different coupling structures. Considering the infinite system-size 
limit, the first step in our analysis has been to reduce the dimensionality of the system. Coarse- 
graining of the network has lead us to the mean-field description. Assuming further that the phases 
in each set of oscillators with the same degree obey a Gaussian distribution with time-dependent 
cumulants, we arrived at a lower-dimensional set of integro-differential equations. The order of this 
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Figure 6. (Color online) Numerical simulations of random binary networks of TV = 10* stochastic rotators. 
Shown are time-dependent probability densities of rotators phases, p{(j>,t). In both panels the rotators are in 
the excitable regime with a = 1.05, D = 0.095 and k = 2. (a): Homogeneous binary network with parameters 
a\ = 0.59, a2 = 0.39 and p — 0.05 leading to the variance Var(Q)/VarMax = 0.0076. (b): Inhomogeneous 
binary network, qi = 0.97, 02 = 0.37, p = 0.05 with the variance Var(Q)/VarMax — 0.0684. 



system equals the doubled number of different degrees occuring in the network. The approach has 
been used to perform numerical bifurcation analysis on regular and binary random networks. The 
calculated state diagrams display the parameter regions where the neuronal network shows syn- 
chronous or asynchronous spiking behavior. Transitions between these regions are typically related 
to passages through the Hopf and saddle- node bifurcations. When the noise intensity is increased 
from zero, the saddle-node bifurcations mark the onset of global firing in the ensemble, whereas 
the Hopf bifurcation leads to the disappearance of the collective mode of synchronous oscillations. 
We have been able to show that changing the network structure alone, without variations in the 
noise intensity and/or the excitation threshold, is sufficient for a switch between the different pat- 
terns of spiking behaviors. In fact, adaptations in the coupling structure play a crucial role for the 
information processing in populations of neurons [3]. 

Our results suggest that deleting uniformly connections of all nodes in the network does not 
lead to new characteristics, since the coupling strength is only rescaled effectively. However, if nodes 
differ in their number of connections, an interesting phenomenon is obtained in the excitable regime. 
We have found out that the critical noise intensity, at which the Hopf bifurcation occurs and small- 
scale oscillations are born, can possess a minimum for a moderate network heterogeneity. It was 
induced by the presence of a small fraction of highly connected nodes. This might have meaningful 
implications for the functioning of neuronal networks, since such hub nodes naturally exist [B^^lD]. 
It is well-known that neuronal networks are highly heterogeneous. Our analysis further suggests 
that this property implies existence of large parameter regions to get global firing. If, however, a 
specific task requests a switching between different spiking patterns, reducing the heterogeneity in 
the connectivities becomes beneficial, since only small noise variations are needed. The optimum 
for such a task would then be given by a moderate coupling variance with a small fraction of hub 
nodes. 
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